Linear Models
Linear models are a basic, but powerful, tool of statistics, and I recommend that everyone serious about visualisation learns at least the basics of how to use them
Diamonds
So far our analysis of the diamonds data has been plagued by the powerful relationship between size and price
It makes it very difficult to see the impact of cut, colour and clarity because higher quality diamonds tend to be smaller, and hence cheaper
carat cut color clarity depth table price
Min. :0.2000 Fair : 1610 D: 6775 SI1 :13065 Min. :43.00 Min. :43.00 Min. : 326
1st Qu.:0.4000 Good : 4906 E: 9797 VS2 :12258 1st Qu.:61.00 1st Qu.:56.00 1st Qu.: 950
Median :0.7000 Very Good:12082 F: 9542 SI2 : 9194 Median :61.80 Median :57.00 Median : 2401
Mean :0.7979 Premium :13791 G:11292 VS1 : 8171 Mean :61.75 Mean :57.46 Mean : 3933
3rd Qu.:1.0400 Ideal :21551 H: 8304 VVS2 : 5066 3rd Qu.:62.50 3rd Qu.:59.00 3rd Qu.: 5324
Max. :5.0100 I: 5422 VVS1 : 3655 Max. :79.00 Max. :95.00 Max. :18823
J: 2808 (Other): 2531
x y z
Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 4.710 1st Qu.: 4.720 1st Qu.: 2.910
Median : 5.700 Median : 5.710 Median : 3.530
Mean : 5.731 Mean : 5.735 Mean : 3.539
3rd Qu.: 6.540 3rd Qu.: 6.540 3rd Qu.: 4.040
Max. :10.740 Max. :58.900 Max. :31.800




ggplot(diamonds %>% filter(carat <= 2), aes(carat, price)) +
geom_bin2d() +
geom_smooth(method = "lm", se = FALSE, size = 2, colour = "yellow")
model <- lm(price ~ carat, data = diamonds)
summary(model)
coef(summary(model))
head(diamonds$price)
head(resid(mod))
summary(resid(mod))
ggplot(diamonds %>% mutate(rel_price = resid(mod)), aes(carat, rel_price)) + geom_point()
ggplot(diamonds %>% mutate(rel_price = resid(mod)), aes(carat, rel_price)) + geom_bin2d()
We can use a linear model to remove the effect of size on price. Instead of looking at the raw price, we can look at the relative price: how valuable is this diamond relative to the average diamond of the same size.

Call:
lm(formula = lprice ~ lcarat, data = diamonds2)
Residuals:
Min 1Q Median 3Q Max
-1.97977 -0.25004 -0.01116 0.24384 1.94623
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.208920 0.002109 5789.0 <2e-16 ***
lcarat 1.696590 0.002078 816.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.377 on 52049 degrees of freedom
Multiple R-squared: 0.9275, Adjusted R-squared: 0.9275
F-statistic: 6.663e+05 on 1 and 52049 DF, p-value: < 2.2e-16
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.20892 0.002108989 5788.9930 0
lcarat 1.69659 0.002078441 816.2798 0
Residuals
residuals are the vertical distance between each point and the line of best fit.
A relative price of zero means that the diamond was at the average price; positive means that it’s means that it’s more expensive than expected (based on its size), and negative means that it’s cheaper than expected.
diamonds2 <- diamonds2 %>% mutate(rel_price = resid(mod))
ggplot(diamonds2, aes(carat, rel_price)) +
geom_bin2d()

xgrid <- seq(-2, 1, by = 1/3)
data.frame(logx = xgrid, x = round(2 ^ xgrid, 2))
color_cut <- diamonds2 %>%
group_by(color, cut) %>%
summarise(
price = mean(price),
rel_price = mean(rel_price)
)
color_cut
ggplot(color_cut, aes(color, price)) +
geom_line(aes(group = cut), colour = "grey80") +
geom_point(aes(colour = cut))

ggplot(color_cut, aes(color, rel_price)) +
geom_line(aes(group = cut), colour = "grey80") +
geom_point(aes(colour = cut))

Exercises
What happens if you repeat the above analysis with all diamonds? (Not just all diamonds with two or fewer carats). What does the strange ge- ometry of log(carat) vs relative price represent? What does the diagonal line without any points represent?
ggplot(diamonds, aes(carat, price)) +
geom_bin2d() +
geom_smooth(method = "lm", se = FALSE, size = 2, color = "Yellow")

ldiamonds <- diamonds %>% mutate(lcarat = log(carat), lprice = log(price))
ggplot(ldiamonds, aes(lcarat, lprice)) +
geom_bin2d() +
geom_smooth(method = "lm", se = FALSE, size = 2, color = "Yellow")

model <- lm(lprice ~ lcarat, data = ldiamonds)
coef(summary(model))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.448661 0.001364691 6190.8959 0
lcarat 1.675817 0.001933806 866.5901 0
ldiamonds <- ldiamonds %>% mutate(rel_price = resid(model))
ggplot(ldiamonds, aes(lcarat, rel_price)) + geom_bin2d()

by_carat <- ldiamonds %>% group_by(lcarat, cut, color, clarity) %>%
summarise(price = mean(price),
rel_price = mean(rel_price))
ggplot(by_carat, aes(lcarat, rel_price, color = cut)) +
geom_point()

ggplot(by_carat, aes(lcarat, rel_price, color = color)) +
geom_point()

ggplot(by_carat, aes(lcarat, rel_price, color = clarity)) +
geom_point()

What does the strange ge- ometry of log(carat) vs relative price represent?
Larger diamonds tend to be cheaper than average because they are typically lower quality
ggplot(by_carat, aes(lcarat, rel_price, color = cut)) +
geom_point()

ggplot(by_carat, aes(lcarat, rel_price, color = color)) +
geom_point()

ggplot(by_carat, aes(lcarat, rel_price, color = clarity)) +
geom_point()

by_carat <- ldiamonds %>% group_by(lcarat) %>%
summarise(price = mean(price),
rel_price = mean(rel_price))
ggplot(by_carat, aes(lcarat, rel_price)) + geom_point()

by_carat_cut <- ldiamonds %>% group_by(lcarat, cut) %>%
summarise(price = mean(price),
rel_price = mean(rel_price))
ggplot(by_carat_cut, aes(lcarat, rel_price)) + geom_point() + geom_smooth(method = "lm")

by_carat_color <- ldiamonds %>% group_by(lcarat, color) %>%
summarise(price = mean(price), rel_price = mean(rel_price))
ggplot(by_carat_color, aes(lcarat, rel_price)) + geom_point() + geom_smooth(method = "lm")

by_carat_clarity <- ldiamonds %>% group_by(lcarat, clarity) %>%
summarise(price = mean(price), rel_price = mean(rel_price))
ggplot(by_carat_clarity, aes(lcarat, rel_price)) + geom_point() + geom_smooth(method = "lm")

I made an unsupported assertion that lower-quality diamonds tend to be larger. Support my claim with a plot.
levels(diamonds$cut)
[1] "Fair" "Good" "Very Good" "Premium" "Ideal"
levels(diamonds$clarity)
[1] "I1" "SI2" "SI1" "VS2" "VS1" "VVS2" "VVS1" "IF"
levels(diamonds$color)
[1] "D" "E" "F" "G" "H" "I" "J"
ggplot(diamonds, aes(log(carat), fill = clarity)) + geom_histogram()

ggplot(diamonds, aes(log(carat), fill = color)) + geom_histogram()

ggplot(diamonds, aes(log(carat), fill = cut)) + geom_histogram()

ggplot(diamonds, aes(carat, clarity, fill = clarity)) + geom_bin2d()

ggplot(diamonds, aes(carat, cut, fill = cut)) + geom_bin2d()

ggplot(diamonds, aes(carat, color, fill = color)) + geom_bin2d()

ggplot(diamonds) + geom_smooth(method = 'lm', aes(carat, color))

ggplot(diamonds) + geom_smooth(method = 'lm', aes(carat, clarity))

ggplot(diamonds) + geom_smooth(method = 'lm', aes(carat, cut))

ggplot(diamonds, aes(log(carat), color=cut)) + geom_density()

ggplot(diamonds, aes(log(carat), color=clarity)) + geom_density()

ggplot(diamonds, aes(log(carat), color=color)) + geom_density()

How do depth and table relate to the relative price?
Depth - no linear relationship


ggplot(diamonds, aes(carat, price)) +
geom_bin2d() +
geom_smooth(method = "lm", se = FALSE, size = 2, color = "Yellow")

ldiamonds <- diamonds %>% mutate(lcarat = log(carat), lprice = log(price))
model <- lm(lprice ~ lcarat, data = ldiamonds)
coef(summary(model))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.448661 0.001364691 6190.8959 0
lcarat 1.675817 0.001933806 866.5901 0
ldiamonds <- ldiamonds %>% mutate(rel_price = resid(model))
by_carat_depth <- ldiamonds %>% group_by(lcarat, depth) %>%
summarise(price = mean(price), rel_price = mean(rel_price))
ggplot(by_carat_depth, aes(lcarat, rel_price)) + geom_point() + geom_smooth(method = "lm")

by_carat_table <- ldiamonds %>% group_by(lcarat, table) %>%
summarise(price = mean(price), rel_price = mean(rel_price))
ggplot(by_carat_table, aes(lcarat, rel_price)) + geom_point() + geom_smooth(method = "lm")

LS0tCnRpdGxlOiAiQ2hhcHRlciAxMSIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6IGRlZmF1bHQKICBodG1sX25vdGVib29rOiBkZWZhdWx0Ci0tLQpgYGB7ciBlY2hvPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgoKIyBNb2RlbGluZwoKKiBVc2luZyBtb2RlbHMgYXMgYSB0b29sIHRvIHJlbW92ZSBvYnZpb3VzIHBhdHRlcm5zIGluIHlvdXIgcGxvdHMKKiBNb2RlbHMgY2FuIGJlIGEgcG93ZXJmdWwgdG9vbCBmb3Igc3VtbWFyaXNpbmcgZGF0YSBzbyB0aGF0IHlvdSBnZXQgYSBoaWdoZXIgbGV2ZWwgdmlldy4KCgoKCiMgTGluZWFyIE1vZGVscwoKCkxpbmVhciBtb2RlbHMgYXJlIGEgYmFzaWMsIGJ1dCBwb3dlcmZ1bCwgdG9vbCBvZiBzdGF0aXN0aWNzLCBhbmQgSSByZWNvbW1lbmQgdGhhdCBldmVyeW9uZSBzZXJpb3VzIGFib3V0IHZpc3VhbGlzYXRpb24gbGVhcm5zIGF0IGxlYXN0IHRoZSBiYXNpY3Mgb2YgaG93IHRvIHVzZSB0aGVtCgoKIyMgIERpYW1vbmRzCgpTbyBmYXIgb3VyIGFuYWx5c2lzIG9mIHRoZSBkaWFtb25kcyBkYXRhIGhhcyBiZWVuIHBsYWd1ZWQgYnkgdGhlIHBvd2VyZnVsIHJlbGF0aW9uc2hpcCBiZXR3ZWVuIHNpemUgYW5kIHByaWNlCgpJdCBtYWtlcyBpdCB2ZXJ5IGRpZmZpY3VsdCB0byBzZWUgdGhlIGltcGFjdCBvZiBjdXQsIGNvbG91ciBhbmQgY2xhcml0eSBiZWNhdXNlIGhpZ2hlciBxdWFsaXR5IGRpYW1vbmRzIHRlbmQgdG8gYmUgc21hbGxlciwgYW5kIGhlbmNlIGNoZWFwZXIKCgpgYGB7ciBlY2hvPUZBTFNFfQpzdW1tYXJ5KGRpYW1vbmRzKQpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhjYXJhdCwgcHJpY2UpKSArIAogIGdlb21fcG9pbnQoKSArIAogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIsIHNlID0gRkFMU0UsIHNpemUgPSAyLCBjb2xvdXIgPSAieWVsbG93IikKCgpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhjYXJhdCwgcHJpY2UsIGNvbG9yID0gY3V0KSkgKyAKICBnZW9tX3BvaW50KCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMikKCgpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhjYXJhdCwgcHJpY2UsIGNvbG9yID0gY2xhcml0eSkpICsgCiAgZ2VvbV9wb2ludCgpICsgCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIiwgc2UgPSBGQUxTRSwgc2l6ZSA9IDIpCgoKZ2dwbG90KGRpYW1vbmRzLCBhZXMoY2FyYXQsIHByaWNlLCBjb2xvciA9IGNvbG9yKSkgKyAKICBnZW9tX3BvaW50KCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMikKYGBgCgoKYGBge3J9CmdncGxvdChkaWFtb25kcyAlPiUgZmlsdGVyKGNhcmF0IDw9IDIpLCBhZXMoY2FyYXQsIHByaWNlKSkgKyAKICBnZW9tX2JpbjJkKCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMiwgY29sb3VyID0gInllbGxvdyIpCmBgYAoKYGBge3J9Cm1vZGVsIDwtIGxtKHByaWNlIH4gY2FyYXQsIGRhdGEgPSBkaWFtb25kcykKc3VtbWFyeShtb2RlbCkKY29lZihzdW1tYXJ5KG1vZGVsKSkKaGVhZChkaWFtb25kcyRwcmljZSkKaGVhZChyZXNpZChtb2QpKQpzdW1tYXJ5KHJlc2lkKG1vZCkpCmdncGxvdChkaWFtb25kcyAlPiUgbXV0YXRlKHJlbF9wcmljZSA9IHJlc2lkKG1vZCkpLCBhZXMoY2FyYXQsIHJlbF9wcmljZSkpICsgZ2VvbV9wb2ludCgpCmdncGxvdChkaWFtb25kcyAlPiUgbXV0YXRlKHJlbF9wcmljZSA9IHJlc2lkKG1vZCkpLCBhZXMoY2FyYXQsIHJlbF9wcmljZSkpICsgZ2VvbV9iaW4yZCgpCmBgYAoKCgpXZSBjYW4gdXNlIGEgbGluZWFyIG1vZGVsIHRvIHJlbW92ZSB0aGUgZWZmZWN0IG9mIHNpemUgb24gcHJpY2UuIEluc3RlYWQgb2YgbG9va2luZyBhdCB0aGUgcmF3IHByaWNlLCB3ZSBjYW4gbG9vayBhdCB0aGUgcmVsYXRpdmUgcHJpY2U6IGhvdyB2YWx1YWJsZSBpcyB0aGlzIGRpYW1vbmQgcmVsYXRpdmUgdG8gdGhlIGF2ZXJhZ2UgZGlhbW9uZCBvZiB0aGUgc2FtZSBzaXplLgoKCmBgYHtyIGVjaG89RkFMU0V9CmRpYW1vbmRzMiA8LSBkaWFtb25kcyAlPiUgCiAgZmlsdGVyKGNhcmF0IDw9IDIpICU+JQogIG11dGF0ZSgKICAgIGxjYXJhdCA9IGxvZzIoY2FyYXQpLAogICAgbHByaWNlID0gbG9nMihwcmljZSkKICApCmRpYW1vbmRzMgoKZ2dwbG90KGRpYW1vbmRzMiwgYWVzKGxjYXJhdCwgbHByaWNlKSkgKyAKICBnZW9tX2JpbjJkKCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMiwgY29sb3VyID0gInllbGxvdyIpCmBgYAoKCgoKYGBge3IgZWNobz1GQUxTRX0KbW9kIDwtIGxtKGxwcmljZSB+IGxjYXJhdCwgZGF0YSA9IGRpYW1vbmRzMikKc3VtbWFyeShtb2QpCmNvZWYoc3VtbWFyeShtb2QpKQpgYGAKCiMjIyBSZXNpZHVhbHMKCnJlc2lkdWFscyBhcmUgdGhlIHZlcnRpY2FsIGRpc3RhbmNlIGJldHdlZW4gZWFjaCBwb2ludCBhbmQgdGhlIGxpbmUgb2YgYmVzdCBmaXQuCgoKQSByZWxhdGl2ZSBwcmljZSBvZiB6ZXJvIG1lYW5zIHRoYXQgdGhlIGRpYW1vbmQgd2FzIGF0IHRoZSBhdmVyYWdlIHByaWNlOyBwb3NpdGl2ZSBtZWFucyB0aGF0IGl04oCZcyBtZWFucyB0aGF0IGl04oCZcyBtb3JlIGV4cGVuc2l2ZSB0aGFuIGV4cGVjdGVkIChiYXNlZCBvbiBpdHMgc2l6ZSksIGFuZCBuZWdhdGl2ZSBtZWFucyB0aGF0IGl04oCZcyBjaGVhcGVyIHRoYW4gZXhwZWN0ZWQuCgoKYGBge3J9CmRpYW1vbmRzMiA8LSBkaWFtb25kczIgJT4lIG11dGF0ZShyZWxfcHJpY2UgPSByZXNpZChtb2QpKQpnZ3Bsb3QoZGlhbW9uZHMyLCBhZXMoY2FyYXQsIHJlbF9wcmljZSkpICsgCiAgZ2VvbV9iaW4yZCgpCmBgYApgYGB7cn0KeGdyaWQgPC0gc2VxKC0yLCAxLCBieSA9IDEvMykKZGF0YS5mcmFtZShsb2d4ID0geGdyaWQsIHggPSByb3VuZCgyIF4geGdyaWQsIDIpKQpgYGAKCmBgYHtyfQpjb2xvcl9jdXQgPC0gZGlhbW9uZHMyICU+JSAKICBncm91cF9ieShjb2xvciwgY3V0KSAlPiUKICBzdW1tYXJpc2UoCiAgICBwcmljZSA9IG1lYW4ocHJpY2UpLCAKICAgIHJlbF9wcmljZSA9IG1lYW4ocmVsX3ByaWNlKQogICkKY29sb3JfY3V0CmBgYAoKYGBge3J9CmdncGxvdChjb2xvcl9jdXQsIGFlcyhjb2xvciwgcHJpY2UpKSArIAogIGdlb21fbGluZShhZXMoZ3JvdXAgPSBjdXQpLCBjb2xvdXIgPSAiZ3JleTgwIikgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91ciA9IGN1dCkpCmBgYAoKCmBgYHtyfQpnZ3Bsb3QoY29sb3JfY3V0LCBhZXMoY29sb3IsIHJlbF9wcmljZSkpICsgCiAgZ2VvbV9saW5lKGFlcyhncm91cCA9IGN1dCksIGNvbG91ciA9ICJncmV5ODAiKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyID0gY3V0KSkKYGBgCgoKCiMjIEV4ZXJjaXNlcwoKV2hhdCBoYXBwZW5zIGlmIHlvdSByZXBlYXQgdGhlIGFib3ZlIGFuYWx5c2lzIHdpdGggYWxsIGRpYW1vbmRzPyAoTm90IGp1c3QgYWxsIGRpYW1vbmRzIHdpdGggdHdvIG9yIGZld2VyIGNhcmF0cykuIFdoYXQgZG9lcyB0aGUgc3RyYW5nZSBnZS0gb21ldHJ5IG9mIGxvZyhjYXJhdCkgdnMgcmVsYXRpdmUgcHJpY2UgcmVwcmVzZW50PyBXaGF0IGRvZXMgdGhlIGRpYWdvbmFsIGxpbmUgd2l0aG91dCBhbnkgcG9pbnRzIHJlcHJlc2VudD8KCgpgYGB7cn0KZ2dwbG90KGRpYW1vbmRzLCBhZXMoY2FyYXQsIHByaWNlKSkgKyAKICBnZW9tX2JpbjJkKCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMiwgY29sb3IgPSAiWWVsbG93IikKCmxkaWFtb25kcyA8LSBkaWFtb25kcyAlPiUgbXV0YXRlKGxjYXJhdCA9IGxvZyhjYXJhdCksIGxwcmljZSA9IGxvZyhwcmljZSkpCgpnZ3Bsb3QobGRpYW1vbmRzLCBhZXMobGNhcmF0LCBscHJpY2UpKSArIAogIGdlb21fYmluMmQoKSArIAogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIsIHNlID0gRkFMU0UsIHNpemUgPSAyLCBjb2xvciA9ICJZZWxsb3ciKQoKbW9kZWwgPC0gbG0obHByaWNlIH4gbGNhcmF0LCBkYXRhID0gbGRpYW1vbmRzKQpjb2VmKHN1bW1hcnkobW9kZWwpKQoKbGRpYW1vbmRzIDwtIGxkaWFtb25kcyAlPiUgbXV0YXRlKHJlbF9wcmljZSA9IHJlc2lkKG1vZGVsKSkKZ2dwbG90KGxkaWFtb25kcywgYWVzKGxjYXJhdCwgcmVsX3ByaWNlKSkgKyBnZW9tX2JpbjJkKCkKCmJ5X2NhcmF0IDwtIGxkaWFtb25kcyAlPiUgZ3JvdXBfYnkobGNhcmF0LCBjdXQsIGNvbG9yLCBjbGFyaXR5KSAlPiUKICBzdW1tYXJpc2UocHJpY2UgPSBtZWFuKHByaWNlKSwKICAgICAgICAgICAgcmVsX3ByaWNlID0gbWVhbihyZWxfcHJpY2UpKQpgYGAKCiMjIFdoYXQgZG9lcyB0aGUgc3RyYW5nZSBnZS0gb21ldHJ5IG9mIGxvZyhjYXJhdCkgdnMgcmVsYXRpdmUgcHJpY2UgcmVwcmVzZW50PwoKTGFyZ2VyIGRpYW1vbmRzIHRlbmQgdG8gYmUgY2hlYXBlciB0aGFuIGF2ZXJhZ2UgYmVjYXVzZSB0aGV5IGFyZSB0eXBpY2FsbHkgbG93ZXIgcXVhbGl0eQoKCmBgYHtyfQpieV9jYXJhdCA8LSBsZGlhbW9uZHMgJT4lIGdyb3VwX2J5KGxjYXJhdCwgY3V0LCBjb2xvciwgY2xhcml0eSkgJT4lCiAgc3VtbWFyaXNlKHByaWNlID0gbWVhbihwcmljZSksCiAgICAgICAgICAgIHJlbF9wcmljZSA9IG1lYW4ocmVsX3ByaWNlKSkKCmdncGxvdChieV9jYXJhdCwgYWVzKGxjYXJhdCwgcmVsX3ByaWNlLCBjb2xvciA9IGN1dCkpICsgCiAgZ2VvbV9wb2ludCgpCgpnZ3Bsb3QoYnlfY2FyYXQsIGFlcyhsY2FyYXQsIHJlbF9wcmljZSwgY29sb3IgPSBjb2xvcikpICsgCiAgZ2VvbV9wb2ludCgpCgpnZ3Bsb3QoYnlfY2FyYXQsIGFlcyhsY2FyYXQsIHJlbF9wcmljZSwgY29sb3IgPSBjbGFyaXR5KSkgKyAKICBnZW9tX3BvaW50KCkKYGBgCgpgYGB7cn0KYnlfY2FyYXQgPC0gbGRpYW1vbmRzICU+JSBncm91cF9ieShsY2FyYXQpICU+JQogIHN1bW1hcmlzZShwcmljZSA9IG1lYW4ocHJpY2UpLAogICAgICAgICAgICByZWxfcHJpY2UgPSBtZWFuKHJlbF9wcmljZSkpCgpnZ3Bsb3QoYnlfY2FyYXQsIGFlcyhsY2FyYXQsIHJlbF9wcmljZSkpICsgZ2VvbV9wb2ludCgpCgoKYnlfY2FyYXRfY3V0IDwtIGxkaWFtb25kcyAlPiUgZ3JvdXBfYnkobGNhcmF0LCBjdXQpICU+JQogIHN1bW1hcmlzZShwcmljZSA9IG1lYW4ocHJpY2UpLAogICAgICAgICAgICByZWxfcHJpY2UgPSBtZWFuKHJlbF9wcmljZSkpCgpnZ3Bsb3QoYnlfY2FyYXRfY3V0LCBhZXMobGNhcmF0LCByZWxfcHJpY2UpKSArIGdlb21fcG9pbnQoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIpCgoKYnlfY2FyYXRfY29sb3IgPC0gbGRpYW1vbmRzICU+JSBncm91cF9ieShsY2FyYXQsIGNvbG9yKSAlPiUKICBzdW1tYXJpc2UocHJpY2UgPSBtZWFuKHByaWNlKSwgcmVsX3ByaWNlID0gbWVhbihyZWxfcHJpY2UpKQoKZ2dwbG90KGJ5X2NhcmF0X2NvbG9yLCBhZXMobGNhcmF0LCByZWxfcHJpY2UpKSArIGdlb21fcG9pbnQoKSAgKyBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKQoKCmJ5X2NhcmF0X2NsYXJpdHkgPC0gbGRpYW1vbmRzICU+JSBncm91cF9ieShsY2FyYXQsIGNsYXJpdHkpICU+JQogIHN1bW1hcmlzZShwcmljZSA9IG1lYW4ocHJpY2UpLCByZWxfcHJpY2UgPSBtZWFuKHJlbF9wcmljZSkpCgpnZ3Bsb3QoYnlfY2FyYXRfY2xhcml0eSwgYWVzKGxjYXJhdCwgcmVsX3ByaWNlKSkgKyBnZW9tX3BvaW50KCkgICsgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikKCgpgYGAKCiMjIEkgbWFkZSBhbiB1bnN1cHBvcnRlZCBhc3NlcnRpb24gdGhhdCBsb3dlci1xdWFsaXR5IGRpYW1vbmRzIHRlbmQgdG8gYmUgbGFyZ2VyLiBTdXBwb3J0IG15IGNsYWltIHdpdGggYSBwbG90LgoKCmBgYHtyfQpsZXZlbHMoZGlhbW9uZHMkY3V0KQpsZXZlbHMoZGlhbW9uZHMkY2xhcml0eSkKbGV2ZWxzKGRpYW1vbmRzJGNvbG9yKQoKZ2dwbG90KGRpYW1vbmRzLCBhZXMobG9nKGNhcmF0KSwgZmlsbCA9IGNsYXJpdHkpKSArIGdlb21faGlzdG9ncmFtKCkKZ2dwbG90KGRpYW1vbmRzLCBhZXMobG9nKGNhcmF0KSwgZmlsbCA9IGNvbG9yKSkgKyBnZW9tX2hpc3RvZ3JhbSgpCmdncGxvdChkaWFtb25kcywgYWVzKGxvZyhjYXJhdCksIGZpbGwgPSBjdXQpKSArIGdlb21faGlzdG9ncmFtKCkKCmdncGxvdChkaWFtb25kcywgYWVzKGNhcmF0LCBjbGFyaXR5LCBmaWxsID0gY2xhcml0eSkpICsgZ2VvbV9iaW4yZCgpCmdncGxvdChkaWFtb25kcywgYWVzKGNhcmF0LCBjdXQsIGZpbGwgPSBjdXQpKSArIGdlb21fYmluMmQoKQpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhjYXJhdCwgY29sb3IsIGZpbGwgPSBjb2xvcikpICsgZ2VvbV9iaW4yZCgpCgoKZ2dwbG90KGRpYW1vbmRzKSArICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBhZXMoY2FyYXQsIGNvbG9yKSkgCmdncGxvdChkaWFtb25kcykgKyAgZ2VvbV9zbW9vdGgobWV0aG9kID0gJ2xtJywgYWVzKGNhcmF0LCBjbGFyaXR5KSkKZ2dwbG90KGRpYW1vbmRzKSArICBnZW9tX3Ntb290aChtZXRob2QgPSAnbG0nLCBhZXMoY2FyYXQsIGN1dCkpCgoKZ2dwbG90KGRpYW1vbmRzLCBhZXMobG9nKGNhcmF0KSwgY29sb3I9Y3V0KSkgKyBnZW9tX2RlbnNpdHkoKQpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhsb2coY2FyYXQpLCBjb2xvcj1jbGFyaXR5KSkgKyBnZW9tX2RlbnNpdHkoKQpnZ3Bsb3QoZGlhbW9uZHMsIGFlcyhsb2coY2FyYXQpLCBjb2xvcj1jb2xvcikpICsgZ2VvbV9kZW5zaXR5KCkKCgpgYGAKSG93IGRvIGRlcHRoIGFuZCB0YWJsZSByZWxhdGUgdG8gdGhlIHJlbGF0aXZlIHByaWNlPwoKCiMgRGVwdGggLSBubyBsaW5lYXIgcmVsYXRpb25zaGlwCgpgYGB7cn0KZ2dwbG90KGRpYW1vbmRzLCBhZXMoZGVwdGgsIHByaWNlKSkgKyAKICBnZW9tX2JpbjJkKCkgKyAKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBzZSA9IEZBTFNFLCBzaXplID0gMikKCmdncGxvdChkaWFtb25kcywgYWVzKGxvZyhkZXB0aCksIGxvZyhwcmljZSkpKSArIAogIGdlb21fYmluMmQoKSArIAogIGdlb21fc21vb3RoKG1ldGhvZCA9ICJsbSIsIHNlID0gRkFMU0UsIHNpemUgPSAyKQoKbW9kZWwgPC0gbG0obG9nKHByaWNlKSB+IGRlcHRoLCBkYXRhID0gZGlhbW9uZHMpCmNvZWYoc3VtbWFyeShtb2RlbCkpCgptX2RpYW1vbmRzIDwtIGRpYW1vbmRzICU+JSBtdXRhdGUocmVsX3ByaWNlID0gcmVzaWQobW9kZWwpKQpnZ3Bsb3QobV9kaWFtb25kcywgYWVzKGRlcHRoLCByZWxfcHJpY2UpKSArIGdlb21fYmluMmQoKQoKCmBgYAoKYGBge3J9CmdncGxvdChkaWFtb25kcywgYWVzKGNhcmF0LCBwcmljZSkpICsgCiAgZ2VvbV9iaW4yZCgpICsgCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIiwgc2UgPSBGQUxTRSwgc2l6ZSA9IDIsIGNvbG9yID0gIlllbGxvdyIpCgpsZGlhbW9uZHMgPC0gZGlhbW9uZHMgJT4lIG11dGF0ZShsY2FyYXQgPSBsb2coY2FyYXQpLCBscHJpY2UgPSBsb2cocHJpY2UpKQoKbW9kZWwgPC0gbG0obHByaWNlIH4gbGNhcmF0LCBkYXRhID0gbGRpYW1vbmRzKQpjb2VmKHN1bW1hcnkobW9kZWwpKQoKbGRpYW1vbmRzIDwtIGxkaWFtb25kcyAlPiUgbXV0YXRlKHJlbF9wcmljZSA9IHJlc2lkKG1vZGVsKSkKCmJ5X2NhcmF0X2RlcHRoIDwtIGxkaWFtb25kcyAlPiUgZ3JvdXBfYnkobGNhcmF0LCBkZXB0aCkgJT4lCiAgc3VtbWFyaXNlKHByaWNlID0gbWVhbihwcmljZSksIHJlbF9wcmljZSA9IG1lYW4ocmVsX3ByaWNlKSkKCmdncGxvdChieV9jYXJhdF9kZXB0aCwgYWVzKGxjYXJhdCwgcmVsX3ByaWNlKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKQoKCmJ5X2NhcmF0X3RhYmxlIDwtIGxkaWFtb25kcyAlPiUgZ3JvdXBfYnkobGNhcmF0LCB0YWJsZSkgJT4lCiAgc3VtbWFyaXNlKHByaWNlID0gbWVhbihwcmljZSksIHJlbF9wcmljZSA9IG1lYW4ocmVsX3ByaWNlKSkKCmdncGxvdChieV9jYXJhdF90YWJsZSwgYWVzKGxjYXJhdCwgcmVsX3ByaWNlKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKQpgYGAKCg==